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We investigate a novel stochastic technique for the global optimization of complex potential energy 
surfaces (PES) that avoids the freezing problem of simulated annealing by allowing the dynamical 
process to tunnel energetically inaccessible regions of the PES by way of a dynamically adjusted 
nonlinear transformation of the original PES. We demonstrate the success of this approach, which 
is characterized by a single adjustable parameter, for three generic hard minimization problems. 
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The development of methods that efficiently determine 
the global minima of complex and rugged energy land- 
scapes remains a challenging problem with applications 
in many scientific and technological areas. In particu- 
lar for NP-hard jjJD problems stochastic methods offer 
an acceptable compromise between the reliability of the 
method and its computational cost. Branch-and-bound 
techniques |j| offer stringent error estimates but scale 
exponentially in their computational effort. In many 
stochastic approaches the computational cost to deter- 
mine the global minimum with a given probability grows 
only as a power-law with the number of variables Q . 

In such techniques the global minimization is per- 
formed through the simulation of a dynamical process 
for a "particle" on the multi-dimensional potential en- 
ergy surface. Widely used is the simulated annealing 
(SA) technique m where the PES is explored in a se- 
ries of Monte-Carlo simulations at successively decreasing 
temperatures. Its success depends often strongly on the 
choice of the cooling schedule, yet even the simplest geo- 
metric cooling schedule is characterized by three param- 
eters (starting temperature, cooling rate and number of 
cooling steps) that must be optimized to obtain adequate 
results. For many difficult problems with rugged energy 
landscapes SA suffers from the notorious "freezing" prob- 
lem, because the escape rate from local minima diverges 
with decreasing temperature. To ameliorate this problem 
many variants of the original algorithm 1^-^ have been 
proposed. Unfortunately these proposals often increase 
the number of parameters even further, which compli- 
cates their application for practical problems. 

In this letter we investigate the stochastic tunneling 
method, a generic physically motivated generalization of 
SA. This approach circumvents the "freezing" problem, 
while reducing the number of problem dependent pa- 
rameters to one. In this investigation we demonstrate 
the success of this approach for three hard minimization 
problems: the Coulomb spin-glass (CSG), the traveling 
salesman problem (TSP) and the determination of low 
autocorrelation binary sequences(LABS) in comparison 



with other techniques. 

Method: The freezing problem in stochastic minimiza- 
tion methods arises when the energy difference between 
"adjacent" local minima on the PES is much smaller 
than the energy of intervening transition states separat- 
ing them. As an example consider the dynamics on the 
model potential in Figure (|f|) (a) . At high temperatures a 
particle can still cross the barriers, but not differentiate 
between the wells. As the temperature drops, the particle 
will get eventually trapped with almost equal probability 
in any of the wells, failing to resolve the energy difference 
between them. The physical idea behind the stochastic 
tunneling method is to allow the particle to "tunnel" |T(]] 
forbidden regions of the PES, once it has been determined 
that they are irrelevant for the low-energy properties of 
the problem. This can be accomplished by applying a 
non-linear transformation to the PES: 



/stun (a) = 1 - exp [~7(/(x) - f )} 



(1) 



where /o is the lowest minimum encountered by the dy- 
namical process so far (see Figure (@)(b) + (<0) §• The 
effective potential preserves the locations of all minima, 
but maps the entire energy space from /o to the max- 
imum of the potential onto the interval [0,1]. At a 
given finite temperature of 0(1), the dynamical process 
can therefore pass through energy barriers of arbitrary 
height, while the low energy-region is still resolved well. 
The degree of steepness of the cutoff of the high-energy 
regions is controlled by the tunneling parameter 7 > 0. 
Continously adjusting the reference energy fo to the best 
energy found so far, successively eliminates irrelevant fea- 
tures of the PES that would trap the dynamical process. 

To illustrate the physical content of the transformation 
we consider a Monte-Carlo (MC) process at some fixed 
inverse temperature (3 on the STUN-PES. A MC-step 
from x\ to X2 with A = /(a^) — f( x i) is accepted with 
probability wi^ 2 = exp [-/3(/stun(^2) - /stun(»i))]- 
In the limit 7A <C 1 this reduces to u)i_>2 ~ exp(— /3A) 
with an effective, energy dependent temperature (3 = 
P^ e i(fo-f(xi)) < rp ne dynamical process on the 
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STUN potential energy surface with fixed temperature 
can thus be interpreted as an MC process with an energy 
dependent temperature on the original PES. In the lat- 
ter process the temperature rises rapidly when the local 
energy is larger than / and the particle diffuses (or tun- 
nels) freely through potential barriers of arbitrary height. 
As better and better minima are found, ever larger por- 
tions of the high-energy part of the PES are flattened 
out. In analogy to the SA approach this behavior can 
be interpreted as a self-adjusting cooling schedule that is 
optimized as the simulation proceeds. 

Since the transformation in equation (|l|) is bounded, 
it is possible to further simplify the method: On the 
fixed energy-scale of the effective potential one can dis- 
tinguish between phases corresponding to a local search 
and "tunneling" phases by comparing /stun with some 
fixed, problem independent predefined threshold ft (see 
Fig. 1(c)). For the success of the method it is essential 
that the minimization process spends some time tunne- 
ls) 4.0 r 



-1.0 - 




FIG. 1. (a) Schematic one-dimensional PES and (b) STUN 
effective potential, where the minimum indicated by the arrow 
is the best minimum found so far. All wells that lie above the 
best minimum found are suppressed. If the dynamical process 
can escape the well around the current ground-state estimate 
it will not be trapped by local minima that are higher in en- 
ergy. Wells with deeper minima are preserved and enhanced, 
(c) After the next minimum to the right has been located, 
wells that were still pronounced in (b) are also suppressed, 
now only the wells around the improved ground-state esti- 
mate and the true ground state are pronounced. Once the 
true ground state has been found (not shown) all other wells 
have been suppressed and will no longer trap the dynamical 
process. The dotted line in (c) illustrates an energy threshold 
< ft < 1 to classify the nature of the dynamics. In order to 
conduct a successful search the dynamical process must ex- 
plore both paths confined to the vicinity of the present well 
(/stun < ft) and paths that escape the well by tunneling the 
barrier (/stun > ft)- Adjusting the temperature to maintain 
a particular average effective energy balances the tunneling 
and the local-search phases of the algorithm. 
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FIG. 2. Average estimated error for the ground-state 
estimates of the Coulomb glass using SA (circles), 
STUN (squares), ST(triangles) and PT(diamonds) for N=100 
(full lines in lower part) and N=500 (dashed lines in upper 
part). 

ling and some time searching at any stage of the mini- 
mization process. We therefore adjust the parameter (3 
accordingly during the simulation: If a short-time moving 
average of /stun exceeds the threshold /thresh ~ 0.03, (3 
is reduced by some fixed factor, otherwise it is increased. 
Following this prescription the method is characterized 
by the single problem-dependent parameter (7). 

Applications: In order to test the performance of this 
algorithm we have investigated three families of compli- 
cated NP-hard minimization problems. For each problem 
we have determined either the exact ground-state energy 
or a good estimate thereof. We computed the average er- 
ror of the various optimization methods as a function of 
the computational effort to determine the computational 
effort required to reach a prescribed accuracy. For the 
applications presented here we have fixed the functional 
form of the transformation and the "cooling schedule" 
for (3 in order to demonstrate that these choices are suffi- 
cient to obtain adequate results. Obviously this does not 
guarantee that these choices are optimal. 

(CSG) The determination of low-energy configurations 
of glassy PES is a notoriously difficult problem. We 
have verified by direct comparison that the method con- 
verges quickly to the exact ground states [jllj for two- 
dimensional short-range Ising spin-glasses of linear di- 
mension 10 to 30 with either discrete or Gaussian distri- 
butions of the coupling parameters. Next we turned to 
the more demanding problem of the Coulomb spin-glass, 
where classical charges {si} with Si — ±1 are placed on 
fixed randomly chosen locations within the unit cube. 
The energy of the system 
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TABLE I. Estimates for the optimal path-length for the traveling salesman problem with N = 20, 50 and 100 sites using 
either only local (left side) or global (right side) moves as described in the text. For global moves both SA and STUN are 
equally efficient to obtain low-energy paths. Using only local moves the existence of barriers hampers the progress of SA. As 
a result SA becomes less efficient than STUN. By virtue of its temperature exchange mechanism PT also allows the random 
walk to cross the barriers, but is less efficient than STUN. The effort is given in thousands of steps, note that the evaluation of 
a local move is much less costly than that of a global move. The path-length indicate the average optimal energy for 20 runs 
and the best energy found. 







Local Moves 


Global Moves | 


N 


Effort 


SA 


PT 


STUN 


SA 


STUN 


20 


50 


4.85 / 3.55 


4.35 / 3.55 


3.60 / 3.55 


3.94 / 3.61 


3.55 / 3.55 


20 


100 


4.52 / 3.58 


4.02 / 3.55 


3.62 / 3.55 


3.93 / 3.55 




20 


500 


4.08 / 3.55 


3.57 / 3.55 


3.55 / 3.55 


3.82 / 3.56 




20 


1000 


4.08 / 3.55 


3.55 / 3.55 








20 


5000 


3.75 / 3.55 










50 


100 


12.5 / 10.61 


13.72 / 12.58 


11.06 / 9.39 


5.74 / 5.65 


5.72 / 5.65 


50 


500 


11.0 / 8.68 


11.55 / 10.65 


8.32 / 5.83 


5.70 / 5.65 


5.67 / 5.65 


50 


1000 


11.0 / 8.84 


10.70 / 9.82 


7.75 / 5.78 


5.68 / 5.65 


5.67 / 5.65 


50 


5000 


9.84 / 8.10 


8.99 / 7.89 


7.16 / 5.78 


5.66 / 5.65 


5.65 / 5.65 


50 


10000 


9.87 / 8.31 




6.70 / 5.72 


5.66 / 5.65 




100 


100 








8.42 / 8.11 


8.40 / 8.01 


100 


500 








8.18 / 8.01 


8.18 / 7.97 


100 


1000 








8.08 / 7.94 


8.03 / 7.95 


100 


5000 








8.01 / 7.94 


8.01 / 7.96 



is minimized as a function of the distribution of the {si}. 

The results of grand-canonical simulations for ten 
replicas of N — 100 and N — 500 charges are shown 
in Figure 2. We first conducted twenty very long STUN 
runs for each replica to determine upper bounds for the 
true ground-state energy. For the same charge distribu- 
tions we then averaged the error of STUN, SA, paral- 
lel tempering (PT) || and simulated tempering(ST) |t]] 
for twenty runs per replica as function of the numerical 
effort. We found that the average STUN energy con- 
verged in 10 6 MC-steps to within 1% of the estimated 
true ground-state energy. Over two decades of the nu- 
merical effort we found a consistent significant advantage 
of the proposed method over the SA approach. Fitting 
the curves in the figure with a power-law dependence we 
estimate that STUN is two orders of magnitude more 
efficient than SA. 

We found no consistent ranking of ST and PT relative 
to SA for the two system sizes considered. Both methods 
offer alternative routes to overcome the freezing problem 
in SA. In PT the configurations of concurrent simulations 
at a variety of temperatures are occasionally exchanged. 
In ST only a single simulation is undertaken, but its tem- 
perature is considered to be a dynamical variable. Tem- 
perature and configuration are distributed according to: 

p(s,T) = e -E(s)/T-g(T) and the weights are Qp _ 

timized for a discretized temperature distribution, such 
that all temperatures are visited with equal probability. 
In both methods, a configuration can escape a local min- 
imum when the instantaneous temperature is increased. 
The choice of the temperature set (along with values for 
g(T)) is system dependent and must be optimized much 



like the annealing schedule in SA. In accordance with 
other studies our results indicate that ST performs signif- 
icantly better than SA for long simulation times. PT was 
successful only for the larger system (N=500), where it 
reached the same accuracy as STUN for 10 6 steps. STUN 
converged faster than any of the competing methods, but 
showed a tendency to level off at high accuracy. In the 
limit of large computational its accuracy was matched by 
ST for N=100 and PT for N=500. 

(TSP) The traveling salesman problem is another ubiq- 
uitous NP-hard minimization problem |l2] , |l4| ]. We have 
investigated the problem in its simplest incarnation: i.e. 
as a minimization of the euclidian distance along a closed 
path of N cities. Using long-range updates, i.e. the rever- 
sal and exchange of paths of arbitrarily length, we found 
that both SA and STUN perform about equally well and 
reach the global optimum for N = 20, 50 and 100 very 
quickly (see right side of Table (||)). 

Nevertheless it is instructive to analyze this model 
somewhat further as it provides insight into the inter- 
play of move-construction and complexity of the min- 
imization problem. The unconstrained TSP is a rare 
instance among NP-hard minimization problems, where 
it is possible to construct efficient "long-range" hops on 
the PES. In most practical applications of minimization 
problems related to the TSP, the construction of global 
moves is severely complicated by the existence of "hard 
constraints" on the routes taken. For such problems, as 
well as the other examples reported here, the alteration 
of just a few variables of the configurations leads to un- 
acceptably high energies in almost all cases. As a result, 
the construction of global moves is not an efficient way 
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to facilitate the escape from local minima. When only 
local moves, i.e. transpositions of two adjacent cities, are 
considered high barriers that were circumvented in the 
presence of global moves hamper the progress of SA. The 
results on the left side of Table (||) demonstrate that in 
this scenario SA performs significantly worse than STUN. 

(LABS) Finally we turn to the construction of low- 
autocorrelation binary sequences HQ. The model can 
be cast as a ground-state problem for a one-dimensional 
classical spin-1/2 chain with long-range four-point inter- 
actions 
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(3) 



and is one of the hardest discrete minimization problems 
known | jj"5| . Even highly sophisticated and specialized 
optimization algorithms [|l4| have failed to find config- 
urations anywhere near (within 20%) the ground-state 
energy that can be extrapolated from exact enumeration 
studies for small systems (N < 50) ]l6| , |l7| ]. The reason 
for this difficulty has been attributed to the "golf-course" 
character of the energy landscape and there is convincing 
evidence that SA will fail to converge to the ground-state 
energy even in the limit of adiabatic cooling Jl3| . The sit- 
uation is significantly improved if the original potential 
energy surface is replaced by a piecewise constant energy 
surface that is obtained by a local minimization of the 
original PES at each point . Obviously the latter sur- 
face preserves all ground-state configurations and ener- 
gies of the original PES, but eliminates many "plateaus" 
of the "golf-course" landscape. Using the modified energy 
surface we are able to compare SA to STUN, since SA 
can now determine the ground state energy of medium 



size systems (N=49) with a large, but finite computa- 
tional effort. Table II summarizes the results for the aver- 
age error of 20 SA and STUN runs for system sizes N=49 
and N=101 as a function of the computational effort. In 
direct comparison we find that STUN is two orders of 
magnitude more efficient than SA. Both methods are at 
least a dozen orders of magnitude more efficient than SA 
on the original PES. 

Discussion: Using three NP-hard minimization prob- 
lems with high-barriers separating local minima we have 
demonstrated that the stochastic tunneling approach of- 
fers a reliable, generic and efficient route for the determi- 
nation of low-energy configurations. One chief advantage 
of the method lies in the fact that only a single param- 
eter must be adjusted to adapt it for a specific problem. 
One of the drawbacks of STUN is that in contrast to e.g. 
PT, no thermodynamic expectation values for the system 
can be obtained from the simulation. Secondly, because 
the non-linear transformation will map any unbounded 
PES onto an interval bounded from above, the dynami- 
cal process in STUN will experience "tunneling" phases 
at any finite temperature. For PES that do not contain 
high barriers, or in the presence of efficient global moves 
that circumvent such barriers, STUN may therefore be 
less efficient than competing methods. In many realistic 
optimization problems where the construction of global 
moves is exceedingly difficult or very expensive the tun- 
neling approach can ameliorate the difficulties associated 
with the existence of high energy barriers that separate 
local minima of the PES. 

We gratefully acknowledge stimulating discussions 
with C. Gros and U. Hansmann. 



TABLE II. Average and best ground state estimates for 
LABS for the N = 49 and N = 101 using SA and STUN 
on the locally minimized PES described in the text. SA now 
systematically approaches the estimated ground-state energy, 
but STUN is about two orders of magnitude more efficient. 
The effort is given in thousands of steps, each step consists of 
a multi-spin flip followed by a local minimization. 
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